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Abstract 

Position determination in biological systems is often achieved through 
protein concentration gradients. Measuring the local concentration of such 
a protein with a spatially- varying distribution allows the measurement of 
position within the system. In order for these systems to work effectively, 
position determination must be robust to noise. Here, we calculate funda- 
mental limits to the precision of position determination by concentration 
gradients due to unavoidable biochemical noise perturbing the gradients. 
We focus on gradient proteins with first order reaction kinetics. Systems 
of this type have been experimentally characterised in both developmen- 
tal and cell biology settings. For a single gradient we show that, through 
time-averaging, great precision can potentially be achieved even with very 
low protein copy numbers. As a second example, we investigate the abil- 
ity of a system with oppositely directed gradients to find its centre. With 
this mechanism, positional precision close to the centre improves more 
slowly with increasing averaging time, and so longer averaging times or 
higher copy numbers are required for high precision. For both single and 
double gradients, we demonstrate the existence of optimal length scales 
for the gradients, where precision is maximized, as well as analyzing how 
precision depends on the size of the concentration measuring apparatus. 
Our results provide fundamental constraints on the positional precision 
supplied by concentration gradients in various contexts, including both in 
developmental biology and also within a single cell. 
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Summary 



In many biological systems gradients of protein concentration provide precise po- 
sitional information. Above a critical concentration a signal can be switched on 
by the gradient, whereas below the threshold the signal is switched off, thereby 
providing position dependent signalling. Such concentration gradients are, how- 
ever, subject to unavoidable noise arising from intrinsic biochemical fluctuations. 
We therefore investigate how precisely a noisy concentration gradient can specify 
positional information. We analyse noisy one and two gradient models, where 
we find that time-averaging of concentration measurements potentially allows 
great precision to be achieved even with remarkably low protein copy numbers. 
We also find that a particular choice of the gradient decay length optimizes 
positional precision. Furthermore, in two-dimensional gradients, such as on a 
cell membrane, we find that positional precision is substantially independent of 
the size of the concentration measuring apparatus. We apply our results to a 
number of examples in both cell and developmental biology, including cell di- 
vision positioning in bacteria and yeast, as well as precision gene expression in 
Drosophila. 

Introduction 

To determine position in a biological system, some component within the system 
must have a non-uniform spatial distribution. Often this is achieved through 
the formation of gradients of protein concentration. Typically a gradient forms 
when a protein is manufactured/injected within a small region, and subsequently 
spreads and decays pQ. By measuring the local concentration, position relative 
to the source can be determined. In developmental biology, where such gradients 
are used to control patterns of gene expression, gradient proteins are called 
morphogens. However, intracellular concentration gradients are also thought to 
be important for organisation inside single cells. 

For a gradient mechanism to be biologically viable, position determination 
must be precise and therefore robust to noise. Variability from one copy of the 
system to another (e.g. from cell to cell or embryo to embryo) will certainly 
compromise positional precision. Production and degradation rates can vary, 
for example, due to different copy numbers of transcription factors or proteases. 
The physical size of the system will also vary and this may affect proper position- 
ing. Most previous analyses of morphogen gradients have focused on robustness 
to changes in these extrinsic factors [2-4] between different copies of the system. 
However, there will also be intrinsic noise affecting the gradient within a single 
copy of the system, for example due to the unavoidably noisy nature of the 
biochemical reactions involved. This dissection of the fluctuations into extrinsic 
or intrinsic mirrors that introduced into the analysis of stochastic gene expres- 
sion [5-7]. However, here intrinsic noise alters not only the overall protein copy 
numbers (similar to [5]), but crucially also the spatiotemporal protein distribu- 
tion. Even if all extrinsic variation could be eliminated, intrinsic biochemical 
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noise would still lead to a fundamental limit to the precision of position deter- 
mination, in a similar way to limits on the precision of protein concentration 
measurement [5J [S]. In this paper, we therefore address the question of how 
precisely a concentration gradient can specify positional information, and cal- 
culate the limits on positional precision for a simple, but biologically relevant, 
gradient formation mechanism with first order reaction kinetics. 

Quantitative measurements, for example on the Bicoid-Hunchback system 
in Drosophila [10] . have shown that remarkable positional precision can some- 
times be obtained. For this reason, understanding the fundamental limits to 
the precision of concentration gradients is clearly an important issue in devel- 
opmental biology. Our results will be equally relevant to gradients that form 
within single cells, where protein copy numbers of a few thousand [11-13] will 
lead to large density fluctuations. The properties of intracellular protein gra- 
dients have been studied by Brown and Kholodenko [14]. Recently a number 
of these gradients have been observed experimentally in both prokaryotic and 
eukaryotic systems. The bacterial virulence factor IcsA forms a polar gradi- 
ent on the cell membrane of Shigella flexneri [15] . MipZ in Caulobacter forms 
polar gradients to aid division site selection [11]. In B. subtilis, the MinCD 
complex also forms polar gradients in order to direct division site selection to 
the mid-plane of the cell [THl Q2] . In E. coli the oscillatory dynamics of the Min 
proteins creates a time-averaged gradient that directs cell division placement 
[18-24]. Using mechanisms of this sort, division site placement in bacteria can 
achieve an impressive precision of ±1% of the cell length [25J [26] . Cell division 
in eukaryotic cells is also believed to be regulated by concentration gradients. 
For example, in fission yeast, the protein Pomlp forms a cortical concentration 
gradient emanating from a cell tip, thereby restricting the cell division protein 
Midlp to the cell centre [2"Tl |2"B] . In eukaryotic cells, gradients of the Ran and 
HURP proteins aid the formation of the mitotic spindle by biasing microtubule 
growth towards the chromosomes [29-33]. Gradients may also play a role in 
the localization of Cdc42 activation, thereby permitting a coupling between cell 
shape and protein activation [311 13S] • 

Suppose that a biological system needs to identify a particular position along 
its length, such as the mid-plane to ensure symmetrical cell division. As con- 
crete examples, MipZ and the MinCD complex act by displacing the essential 
cell division protein FtsZ from the cell membrane. Since the concentrations 
of MipZ/MinCD are higher near the cell poles, FtsZ accumulates near the cell 
centre. Below some critical threshold of MinCD or MipZ concentration, enough 
FtsZ will presumably accumulate to form the division apparatus. The locations 
where the concentration gradient crosses these thresholds mark positions within 
the cell. In our analysis we will simply postulate the existence of such well- 
defined critical thresholds, where the gradient sharply switches a downstream 
signal from on to off. Clearly any real gradient cannot act as such a sharp switch 
- in reality a certain amount of smearing is inevitable. Furthermore, there will 
be additional noise in the process of actually measuring the concentration due 
both to the binding of the gradient proteins to the receptor molecules [8] [9] , and 
also to the downstream reactions that process this incoming signal [5-7,36-38]. 
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In general, the noise of the output signal of a processing network can be written 
as the sum of a contribution from the noise in the input signal plus a contribu- 
tion from the reactions that constitute the processing network. We assume here 
that the detector and the processing network are ideal and do not add any noise 
to the gradient input signal. As a result, our calculated variation constitutes 
a lower bound; any real gradient signalling system will inevitably have a lower 
precision. 

We first consider a system with a single planar morphogen source and lin- 
ear degradation, thereby producing an exponentially decaying average concen- 
tration profile. While this model is very simple, it remains biologically rele- 
vant in both developmental and intracellular contexts. Gradients of Bicoid in 
Drosophila and IcsA in Shigella have been quantitatively measured and shown to 
fit this exponential decay profile on average to high accuracy [lOl [T5] . We then 
calculate the expected distribution of positions where a noisy gradient crosses a 
concentration threshold. With typical cellular copy numbers of a few thousand 
proteins, the system will be unable to identify the correct threshold position 
from a single measurement. In order to achieve reliable position determina- 
tion the concentration must be averaged over time. We show that by averaging 
measurements, a biological system is able to achieve precision in position de- 
termination of a few percent of the system size even with hundreds of protein 
copies, a result we verify by computer simulations. Furthermore, we find that 
the precision of position determination is maximised when a particular choice of 
the gradient decay length is made. We also show how the precision depends on 
the detector size (i.e. the volume over which the density measurement is made). 
For a two dimensional gradient (e.g. on a membrane), the precision possible 
after a certain averaging time depends only very weakly on the detector size. 
We relate all these results to experimental measurements of gradients in Shigella 
and fission yeast. 

We also consider the ability of gradients from two poles to identify the centre 
of the system, as in the MipZ and Pomlp gradients discussed above. Related 
designs have also been proposed for the control of hunchback positioning in 
Drosophila [3, 4, 39J. As before, we find that the precision of the system can be 
optimised by a particular choice of the decay length. However, if the threshold 
position is set at the system centre, time-averaging improves precision more 
slowly than in the single-source model. For subcellular gradients we find that a 
few thousand copies of the gradient proteins may therefore be required for high 
precision. Our results strongly constrain the possible concentrations of gradient 
proteins in two gradient systems. 

Results 

Single Gradient Model 

We consider a protein gradient which is used to determine a particular position 
along the length of a cylindrical system. The system will have dimension d = 2 
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if the gradient is restricted to the membrane, or d = 3 if the gradient is in the 
cytoplasm. We choose the x-axis along the long axis of the system. Position in 
the remaining coordinates is denoted by the vector y. For a membrane system, 
periodic boundary conditions are appropriate in the y-direction. Otherwise, 
zero- flux boundaries are used throughout. The system length is L, and the size 
of the system in the remaining directions is taken to be L± (so L±_ = 2irr, where 
r is the system radius, for the 2d membrane case). A source on the x = plane 
produces proteins at rate J per unit area, which then diffuse with diffusion 
constant D, and are degraded uniformly at rate p. Neglecting fluctuations, the 
protein concentration p(x, y, t) will be described by 

^=DV 2 p- pp + J6(x). (1) 

If L 3> A = y/D/fi, the characteristic decay length of the gradient, we find that, 
at steady state, the density is 

p(x) = ^exp(-x/X). (2) 

Symmetry dictates that the average density is independent of y. Gradients 
with the form ([2|) have been found to accurately fit quantitatively measured 
concentration profiles in both developmental [10] and subcellular p~5] systems. 

While we have outlined the model in terms of production and degradation, 
([l} could equally apply to other mechanisms in which the active protein orig- 
inates in a single location, but deactivation occurs uniformly throughout the 
system. The same equation would therefore describe a protein which is phos- 
phorylated by a polar-localised kinase and dephosphorylated by a uniformly 
distributed phosphatase, or a protein which is activated by being injected into 
the membrane at a pole and deactivated when it dissociates. These biochemical 
details do not affect the behaviour of the model. 

We suppose that signalling is active where the local gradient protein concen- 
tration is above some threshold value, pr, and inactive otherwise. The average 
concentration profile for a single gradient, @, suggests that the system will 
be divided into a region < x < xt where signalling is active, and a region 
xt < x < L where signalling is not active, with pt = p(xx)- However, noise in 
the local protein concentration will cause this threshold position to fluctuate. 
This noise may come from intrinsic fluctuations in the diffusion, injection and 
decay processes, or from extrinsic factors which produce systematic changes in 
the boundary position when comparing one copy of the system to another. Here 
we consider only intrinsic biochemical fluctuations. 

Protein production and degradation events are considered to be single molecule 
reactions with a fixed probability per unit time, and hence will be Poisson pro- 
cesses. We also assume that the hopping of proteins in or out of a particular 
region of space is governed by Poisson statistics, thereby generating a diffusive 
process for molecular transport. Since the system is linear, the instantaneous 
fluctuations in molecular number, n, within a volume (Ax) d centred on the 



5 



position (x, y) should also obey Poisson statistics, with 



(n(x) 2 )-(n(x)) 2 = (n(x)). 



(3) 



In terms of protein density, this becomes 



((Ap(x)f) = (p(x) 2 ) 



(4) 



(A X y 



This relation can also be established using more elaborate field theoretic tech- 
niques (see [10]). From this expression for the variation in the density we can 
compute the width of the threshold position distribution by expanding about 
the average threshold position xt- To leading order, this width is given by 



where p'(xt) denotes the first derivative of the density at x = xt- 

Here we identify (Ax) d as the size of the region in which the concentration 
is being measured. For subcellular gradients involved in positional information, 
this volume will be determined by the size of an individual receptor or protein 
with which the gradient protein interacts, an example being the interaction 
between the MinCD and FtsZ proteins in B. subtilis. The size of the detector, 
Aa;, will then be on a molecular scale. This conclusion still holds even if the 
gradient proteins bind cooperatively to the "detection" protein/receptor due to 
the close physical proximity of the bound molecules. In contrast, however, the 
cellular length scale will be much larger, 1/im or bigger. 

Throughout the following analysis we will focus on subcellular gradients. 
However, our model can equally be applied to developmental biology, and we will 
consider these systems further in the Discussion. As concrete examples we first 
consider the IcsA polar gradient on the membrane of the rod-shaped bacterium 
Shigella (L 3/im, L± « 3^tm) |15j . IcsA is exported to the outer membrane 
at a single pole, after which it diffuses and undergoes uniform proteolysis by the 
protease IcsP, thereby forming an exponential gradient exactly as in our model 
[15]. Outer membrane IcsA is then able to recruit actin nucleation factors. 
However, a critical concentration of IcsA is likely needed for actin nucleation: 
in this way a comet-like actin tail is generated at only one cell pole thereby 
generating unidirectional motility of the pathogen. A cell will typically have 
a few thousand copies of IcsA [T2], forming a gradient with A « 0.5/im [T5] . 
We take the detector size to be Ax — 0.01/im, consistent with an interaction 
between IcsA and actin nucleation proteins. For diffusion on the cell membrane, 
we take D = l/im 2 s _1 . On the membrane of a cell of this size, there would be 
approximately LL±/ (Ax) 2 ~ 10 5 potential detector sites, many more than the 
typical copy number. Even near to the source pole, detector sites will typically 
be unoccupied. A detector region at a distance x — 0.5/im from the highly- 
occupied pole will have average occupancy of (n) ~ 10 _1 . In the cytoplasm 
of a similarly sized bacterium, the number of potential detector sites will be 





6 



~ 10 6 , again much larger than the protein copy numbers typically supported by 
bacteria. 

Similar estimates can be made for single polar gradients in fission yeast 
(L = lOjum, L± = 6/im), such as for Pomlp [571 [55]. Here we assume a total 
of 2000 protein copies (this concentration has not yet been measured but this 
number is plausible [IE]). We also take D = l[j,m 2 s~ 1 and a decay length 
of A = 2/im, parameters that are approximately consistent with the Pomlp 
gradient imaged by Padte et al [25]. We again assume that Ax — O.Ol^m 
corresponding to a molecular sized detector, as would be the case if the gradient 
protein interacted with other membrane proteins (such as Midlp) [2Z112S]- The 
typical occupancy of a Ax = 0.01/im site is then (n) ~ 10~ 2 at x — 2/xm from 
the source. 

As we have seen for both fission yeast and Shigella, average detector site 
occupancies that are very much less than one ensure that the threshold concen- 
tration must necessarily be less than one protein per site. Since most regions 
will be devoid of any copies of the protein, a single instantaneous measurement 
of the protein density cannot give a good estimate of the local average con- 
centration. Additionally, multiple positions where the concentration crosses px 
will be observed simultaneously in such a measurement since the concentration 
will be above the threshold everywhere there is a protein molecule present, and 
below the threshold where there is no protein molecule. In order to reliably 
determine the average concentration profile the system must therefore integrate 
the measured concentration over time. 

The noisy concentration profile provided by the gradient protein forms the 
input signal that is then time-averaged by a downstream signal processing net- 
work. In general, the mechanism for time averaging is provided by the lifetimes 
of the states in the processing network. For instance, in the case of gene expres- 
sion, fluctuations in the occupancy of the promoter by a gene regulatory protein 
can be filtered by the lifetime of the mRNA transcript, provided that lifetime 
is much longer than the timescale of fluctuations in the promoter occupancy 
[7J [5]. Similarly, for subcellular gradients, as in Shigella, fluctuations in the 
gradient can be filtered by the lifetime of activated receptors/detector proteins 
or their downstream products. Provided this time scale is much longer than the 
sub-millisecond timescale of the gradient fluctuations, then good time-averaging 
can be achieved. Importantly, the reactions in the downstream network not only 
time-average the noise of the input signal, but also add further noise to the sig- 
nal [5-7,36-38]. Here, we focus exclusively on noise in the concentration gradient 
and do not model the downstream reactions explicitly, but simply assume they 
are noiseless and model them with an effective averaging time. In essence we 
assume that the detector and the network that the process the gradient signal 
are ideal and do not add further noise, and are thus able to time-average the 
gradient signal in the best possible way. Our results thus provide a lower bound 
to the output noise set by the Poissonian fluctuations of the signalling molecules. 

We suppose that averaging over a time-interval r we can take N T = r/rind 
independent measurements of the concentration. In our ideal case, we then 
expect that the fluctuations in the concentration will decrease according to 
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1 / y/N T . Since the width varies linearly with Ap according to , the width will 
also decrease as 

w(t) ~ w ^l^^. (6) 

The time-scale Ti„d on which independent measurements can be made is set in 
our ideal case solely by the reaction-diffusion dynamics of the gradient proteins, 
as discussed in the Appendix. For cellular parameter values, the typical reaction 
timescale, will be much longer than the typical timescale for diffusion 

between detector sites, (Ax) 2 /D. Assuming a molecular sized detector, this 
latter timescale will be of order 10~ 4 s, whereas effective protein lifetimes will 
typically be seconds or longer. The Damkohler number for the system, the ratio 
of the diffusive and reaction timescales, will therefore be Da ~ (Ax) 2 /A 2 ~ 
I0~ 4 . Since Da <C 1, the averaging time-scale is dominated by diffusive motion. 
In d = 3 we find Ti n d ~ (Ax) 2 /D. However, in d = 2, density correlations 
decay away more slowly, leading to the appearance of logarithmic corrections 
that are weakly dependent on the parameters A and Ax. For long averaging 
times, t 3> l/p, the width determined from time-averaged measurements will 
be 



w(t) = k 2 d 
d = 2, and for d = 3 



w(t) = k 3d 



A 



tJ(Ax) 



exp (xt I A) 



1/2 

(7) 



1/2 

(8) 



where k 2 d, k 3 d and a are constants. 

As we have discussed above, Aa; will be set by the concentration detection 
mechanism. However, in a subcellular context, Aa; also sets the highest possi- 
ble resolution of the system. Once w ~ Aa; the cell cannot resolve the target 
position with any higher precision. Equation (0 suggests that in d = 2, preci- 
sion dependends only very weakly on the detector size, through the logarithmic 
correction factor. Reducing the detector size will increase the number of inde- 
pendent measurements made in a given averaging time. However, since fewer 
proteins will be measured by each detector over one averaging period, reducing 
Aa; will therefore increase the instantaneous density fluctuations. In d = 2 these 
two effects will largely cancel. Hence, even if we have over/underestimated the 
detector volume, this will have little effect on the precision of two dimensional 
gradients, such as IcsA in Shigella or Pomlp in fission yeast. In three dimen- 
sions, however, w varies as (Aa;) -1 / 2 . Since increasing Ax reduces w in both 
d = 2 and d = 3, an optimal strategy would be to choose Aa; to match the 
desired precision in order to minimise the required averaging time. 

Intriguingly, from equations {7J and J5J) we find that there exists an optimal 
decay length such that precision is maximised. This result can be understood 
as follows: for fixed xt, and for A ^> xt, the value of the | (p'(xt)} \ tends to 
a constant J/D, independent of xt- However, as A increases, (p(xx)) increases 
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and therefore the absolute size of the fluctuations in the density also increases. 
Therefore, for large and increasing values of A, w oc {\/ p(xr)) / (p'{ x t)) must 
be increasing. Now if A is small (A -C xt) and decreasing, when computing 
the width oc (a/ p(xt)} / {p' (xt)} the presence of the square root means that the 
numerator decreases much more slowly than the denominator. Hence the width 
must again increase as A is decreased for small A. Combining these results for 
small and large A, the width must have a minimum, optimum value as a function 
of A. This occurs at A = xt in d = 3. In d = 2, the optimal decay length is 
given approximately by 



where we have retained the first order logarithmic correction. 

In order to examine the biological impact of equation (|7|) we again consider 
the Pomlp membrane gradient in fission yeast [27l [28] , using the parameters 
described earlier. Simulations of this example system were performed as de- 
scribed in the Methods with on average 100 proteins in the system. Figures UK 
and B show how the measured threshold position, x, and width, w, vary with 
averaging time. For long averaging times the simulation data gives excellent 
agreement with (JTJ), with the constants kid = 0.40 ± 0.02 and a = 2.5 ± 0.8. 
Figure \Vp shows the w ~ r -1 / 2 behaviour predicted in ([7|), and figure QJ) con- 
firms that the width has a minimum as a function of A. The simulation results 
are consistent with the position of the minimum predicted by Figure [If] 
shows that the distribution of measured threshold positions is Gaussian to a 
good approximation. 

Since the averaging timescale Ti n d in a subcellular system is of order ~ 10~ 4 s, 
time-averaging over a period of minutes can achieve great precision even with 
very few copies of the gradient protein. With the parameter values given above, 
equation ([7|) predicts that the position xt — 2/im can be located to within 
±0.5/im within an averaging time r = 60s even if the system contains on average 
only about 20 copies of the protein. ±0.1/^m precision can be achieved in the 
same averaging time with around 400 copies of the protein, a remarkably high 
level of precision for such a low concentration. In vivo Pomlp gradients may be 
formed by a few thousand protein copies, allowing for even greater precision. 

However, we can see in figure [TJ3 that for averaging times of less than about 
a second, the simulation results are not consistent with ([7]). In this regime 
both w and x are equal to A. As discussed above, at very short averaging 
times the presence of a particle at any position will cause the time-averaged 
concentration to be above pt at that point and hence generally will generate 
a threshold crossing. The probability distribution of threshold measurements, 
p{x), will therefore follow the probability distribution of particles. Assuming 
L>Awe have 





p(x)dx — A 1 exp(— x/X)dx. 



(10) 
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The cell will on average estimate the threshold position to be 

x = I xp(x)dx s» A, 
Jo 



(11) 



and measurements will be distributed about this position with variance 



The system is therefore unable to resolve the correct threshold position at these 
short time scales if this is different from A. 

Associated with the average concentration at the threshold is a length scale, 
/ ~ Pj- 1 ^! the typical distance between proteins at this position. The average 
time for a protein to diffuse this distance will scale as l 2 /D. In two dimensions, 
this time is given by 



Since r x is the timescale on which a diffusing particle first arrives at xt, if 
r <C t x there will generally be no particles detected at xt in the averaging 
period. The system therefore cannot reliably estimate the mean concentration 
at xt, and hence cannot precisely identify the threshold position. For averaging 
times much greater than t x , on average at least one particle will be detected 
at xt- The time- averaged concentration profile will then approach and 
x will approach xt- Hence t x determines the cross-over time between the two 
observed regimes of constant w and w oc t~ 1//2 . Figure[lp shows that the scaling 
in equation (TT3")) is also reproduced in our simulations. For the parameter values 
above, r x = 0.3s, and for a more realistic copy number of 1000, t x = 0.03s. 
These timescales are extremely short compared to cell cycle timescales, but do 
nevertheless show that some sort of time averaging is probably essential: a single 
instantaneous measurement is unlikely to provide precise positional information. 
In fact, as we have seen, averaging over much longer times (tens of seconds) may 
be necessary if very high (1%) precision is required. 

Simulations of the model in three dimensions were also performed (data not 
shown). Similar behaviour was observed in this case, and equation JHJl gave 
good agreement with the observed width at long averaging times. 

Oppositely directed gradients 

In order to reliably locate the centre of a system, the mechanism responsible 
must incorporate information about the overall system size so that the iden- 
tified position can scale correctly. A single gradient characterised by a fixed 
decay length cannot achieve this. We therefore examine a system where protein 
gradients are produced by sources at both ends, and where the central position 
is identified as a concentration minimum. 

We modify our earlier model by adding an additional planar source at x = L. 
This addition is appropriate for modelling cell division inhibitors, such as MipZ 




(12) 



((p(x T )) D)- 1 = (JA)- 1 exp(x T /A). 



(13) 



10 



in Caulobacter, that are injected into the membrane near both cell poles. How- 
ever, our model would apply equally if the two sources are of different repressor 
proteins (as may be the case in fission yeast [UJ HB]), although we do assume 
that J, D and p are the same for both gradients. In this scenario, signalling 
activity will be determined by the total concentration. Without fluctuations, 
this will be described by 

^ = DV 2 p- pp + JS(x) + J6(x - L). (14) 

The steady-state solution is now 

JA cosh((x-£/2)/A) 
P[X) ~ ~D sinh(L/2A) ' (15) 

which has the expected minimum at x = L/2. 

We then suppose that the cell compares the concentration to a threshold 
value corresponding to the minimum of the average profile, p m in = p(L/2) = pr- 
Positions where the concentration is at or below the threshold are identified as 
being at the centre of the cell. While the average steady-state density profile 
would never extend below p m in, fluctuations ensure that the concentration in 
the region around the centre spends a significant amount of time at or below the 
threshold. Around point(s) where (p(x)) = pr, noise in the protein concentra- 
tion will lead to a distribution of threshold crossing positions. We consider an 
expansion of the density fluctuations about xt = L/2, giving, to leading order 

Ap(x T ) = ±\(p"(x T ))\w 2 , (16) 

since any first order term proportional to (p') vanishes at xt = L/2. The width 
is therefore given by 

2 2Ap(L/2) 



- - wm- <i7) 



Substituting in fig]) gives 



/ 4JA^sinh(L/2A) \ 1/4 

w °-[ JiKxT ) ■ (8) 

As in the single gradient model, the typical occupancy of the threshold 
region will be much less than one. For example, if we take the parameter 
values considered previously for the Pomlp gradient in fission yeast, with 2000 
protein copies, the average occupancy of a detector site at x = L/2 will be 
(n{L/2)) ~ 10~ 3 . We assume here that Pomlp forms a gradient from both poles. 
In fact it may only form a single gradient with another hitherto unidentified 
protein forming the second polar gradient [57J[2H]- However, as discussed earlier, 
this detail does not affect our calculations. As a second example, MipZ in 
Caulobacter (L = 2.5pm, L± = 2pm) is typically present at about 1000 copies, 
and forms two polar gradients with a decay length A ps 0.25/im The average 
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occupancy at the centre of this system would be approximately (n(L/2)) ~ 10~ 3 . 
Averaging measurements of the concentration over time is therefore required in 
both cases to obtain precise positional information. Since the width now goes 

1/2 

as (Ap) ' , as shown in (TTTj) , we expect 



w(t) 



w 



1/4 



= < 



k 2 d 



A J sinh(L/2A) (In 



(Ax) 2 
i/4 



1/4 



in d = 2 



in d = 3 



(19) 



where fc2d, <5 and £;3d are constants. Averaging proceeds much more slowly than 
previously, with a r -1 / 4 dependence. This follows directly from the vanishing of 
the first derivative at the average threshold position. In d = 3, and for A<L. 
equation (|19p predicts that w will be minimised when A « L/6 is chosen. In 
d = 2 logarithmic corrections again alter this result slightly, with the optimal 
decay length now occurring at 



1 - 



1 



31n(L/6(Ax)) 



(20) 



where we have included the leading logarithmic correction. This result arises for 
similar reasons as in the single gradient model. For the Pomlp gradient imaged 
by Padte et al [2S], the decay length is observed to be 1 — 1.5/im, comparable 
to this optimal decay length of about 1.5/im for a 10/im cell. 

We simulated our model in two dimensions with representative parameter 
values for fission yeast membrane gradients. We used /j, = 0.36s -1 chosen to 
give A = 1.67/im, and J — 6/im~ 1 s~ 1 giving on average 200 protein copies 
in total. Figure shows the results of these simulations. Again we observe 
two distinct regimes. At averaging times longer than about a second, there is 
excellent agreement with equation (| 19[) . as we can see in figure [2]C. Fitting to 
the simulation results we find k 2 d — 0.63 ± 0.02 and a — 2.5 ± 1.0. Figure 
confirms the existence of the optimal decay length in our simulations. 

Since the width decays as r -1 / 4 for this system, longer averaging times 
and/or higher protein copy numbers are required than in the single gradient 
model to achieve high precision. Intrinsic biochemical noise may therefore 
strongly constrain systems of this type. In order for the yeast-membrane gradi- 
ent considered above to achieve precision of ±5% of the cell length after aver- 
aging for one minute, about 800 protein copies are required. Therefore, in the 
absence of any other positioning mechanisms, the Pomlp gradient will require 
~ 1000 protein copies or more to precisely direct the location of cell division. 
We estimate that the MipZ gradient in Caulobacter, with 1000 protein copies, 
would be able to locate the cell centre to within ±5% of L after approximately 
t = 2s. However, since precision only improves as t -1 / 4 , averaging over t = 20 
minutes would be required for the same system to achieve ±1% accuracy. 
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Discussion 



Noise in biochemical processes within a cell will lead to fluctuations in pro- 
tein concentration gradients, and hence also to variation in the position where 
these gradients cross a particular threshold value. These fluctuations therefore 
place a limit on the potential precision of position determination mechanisms 
relying on concentration gradients alone. In subcellular systems with protein 
copy numbers in the thousands, this noise will be sufficiently large that position 
cannot be determined reliably from a single measurement of the density profile. 
In order to determine position to within a few percent, a precision achieved 
by some subcellular systems, the protein concentration must be averaged over 
time. For a single subcellular membrane gradient, we have seen that by aver- 
aging over a period of a minute, excellent precision can potentially be achieved 
with only a few hundred protein copies. This remarkable precision is due to 
the sub-millisecond diffusive time-scale on which time-averaging occurs. Precise 
identification of the cell mid-plane by gradients emanating from both poles re- 
quires longer averaging times or higher copy numbers, since larger fluctuations 
result from the vanishing first derivative of the average concentration at the sys- 
tem centre. Intrinsic biochemical noise may therefore be a strong constraint on 
subcellular two-gradient positioning systems, dictating that the copy numbers 
be sufficiently high to suppress fluctuations. 

So far we have focused almost exclusively on fluctuations in subcellular gra- 
dients, however our results are also applicable to developmental biology and we 
wish to briefly comment on this application. Here the appropriate length scales 
are usually much longer, on the order of hundreds of micrometers in Drosophila. 
Moreover, the gradients affect patterns of gene expression through the binding 
of gradient molecules to DNA regulatory sequences inside individual nuclei. For 
example, in Drosophila, where exponential gradients have been quantitatively 
measured for Bicoid [10] , Bicoid binds cooperatively to hunchback regulatory 
DNA. In this case we again expect molecular-scale effective measuring volumes, 
with Ax ~ 0.01/im being a reasonable order of magnitude. We next assume 
purely Poisson statistics for the fluctuations: this is a stronger assumption than 
for our earlier subcellular gradients, as there will be additional complications 
arising, for example, from the import/export of morphogens from nuclear com- 
partments. However, if diffusive noise is dominant then Poisson statistics will 
be retained and we can expect our earlier analysis to apply, although with one 
important distinction. Instead of Ax setting the maximal possible precision, 
this will now be set by the size of individual nuclei (prior to ccllularization), 
since we expect relatively homogeneous gene expression within a single nuclear 
volume. A single nucleus in Drosophila has a length scale of around 10/im, still 
much smaller than the decay length of the gradient of A ~ 100/im, allowing for 
high precision gene expression [10]. Using the Drosophila Bicoid gradient as an 
example, we use L — 500/im, Lj_ = 100/im, and estimate D — 10/im 2 s~ 1 and 
/i = 10 _3 s~ 1 , giving A = 100/im, consistent with experiment |10j . Assuming a 
high copy number of 10 7 per embryo (we are not aware of experimental con- 
straints on this figure), gives J ~ l/im _2 s _1 . For a single gradient in three 
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dimensions, wc find that about a 5 minute averaging time is required to bring 
the error down to plus or minus a single nuclear length. For a two gradient 
model in three dimensions, longer averaging times on the order of an hour are 
required to reduce the centre-finding positional error to plus or minus about 
2 nuclear lengths. Since gene expression may need to be controlled on shorter 
timescales than this, other designs, for example using interacting gradients pHH], 
may be required for high precision centre finding (see also below). The effects 
of the optimum gradient length scale will also be interesting to probe in a de- 
velopmental biology context. However, our simple analysis may be complicated 
by the multiple roles played by many morphogens: for example, Bicoid not 
only activates hunchback, but it also helps to regulate pair-rule genes, such as 
Even-skipped. Nevertheless, it is interesting to note that the Bicoid gradient 
length scale A ~ 100/im [10] is not too far away from the L/6 optimum for a 
two gradient case, and in a single gradient context will offer maximal precision 
well into the anterior half of the embryo. 

Up to this point we have only considered systems with first order degra- 
dation. Morphogen gradients with nonlinear decay have also been proposed 
[2|. This nonlinearity will lead to non-Poissonian density fluctuations, which 
may significantly change the observed behaviour. England and Cardy [41j have 
previously calculated the response of a gradient with nonlinear decay to one 
source of biochemical noise, namely a fluctuating production rate. However, 
they calculated the change to the average gradient, while fluctuations about 
this average may also be important. It would certainly be of interest to com- 
pare the performance of linear and nonlinear degradation mechanisms in more 
detail. Centre-finding mechanisms with interactions have also been proposed 
[3 |4] . In these models position is determined from the combined gradient of 
two proteins, which will be steep around the system centre due to an interac- 
tion between the two gradients. These mechanisms may therefore be able to 
achieve greater precision for mid-point determination than the noninteracting 
mechanism considered here. 

Throughout this work we have assumed that the gradient protein concentra- 
tion fluctuates about a steady-state profile, and hence averaging over a longer 
time will give a more precise estimate of the average profile. For a subcellular 
system, the steady-state gradient will develop over timescales of less than about 
a minute, due to the micrometer length scales involved. This timescale is short 
compared to the cell cycle time, which ranges from tens of minutes up to many 
hours. For this reason we expect that subcellular gradients will be in steady- 
state and therefore that our analysis will be directly applicable. However, in 
developmental biology, the effective lifetimes will likely be much longer, and 
the gradient may take hours to fully reach steady-state. Moreover, a number 
of developmental biology systems are known to respond to a morphogen gra- 
dient that has not reached steady-state [42-44]. A further complication is the 
possibility of gradient formation by non-Fickian diffusion |45) . where there is 
no steady-state at all. The model considered in this paper does not take into 
account time-varying average gradients. If the average gradient is evolving, a 
longer averaging period will not necessarily lead to improved precision. Clearly, 
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more work will be required to understand how such dynamically evolving sys- 
tems are able to yield precise positional information and filter out fluctuations. 
Nevertheless, we do note that two gradient systems of the kind analyzed here are 
naturally able to locate the system centre even without being in steady-state, 
due to the symmetry of the system [5] . The positional variations in such a non- 
steady-state scenario will not be the same as calculated here, but our analysis 
does form a first step towards the analysis of these more complex systems. 



Methods 
Calculation of r ind - 

We have assumed in our analysis that during the time-averaging process we are 
taking independent measurements at intervals of However, in both real 

biological systems and in our simulations, measurements can generally be taken 
at much shorter intervals than this, leading to correlations between consecutive 
measurements. For a series of correlated measurements taken at time intervals 
St over a period < t < r, with r 3> St, the expected error for the time-averaged 
concentration at position x, (Ap(x, t)) 2 , is given by [46] 



(Ap(x,r)) 2 = -(Ap(x,0)) 2 

T 



(21) 



i+ ir( i -7) c(f) " i 

where (A/j(x, 0)) 2 is the variance of a single measurement, 

(Ap(x,0)) 2 = (p(x,0) 2 )-(p(x,0)) 2 , (22) 

and C(t) is the normalized density correlation function, 

c(t) = (p(x,t)p(x,0))-(p(x,0)) 2 
(p(x,0) 2 )-(p(x,0)) 2 

We therefore define the timescale Tind to be 

T ind {T)=2j (l - t\ C(t)dt, 
and assuming T in( i ^> St we recover 



(24) 



Ap(x,T) = Ap(x,0)(^^) 1/2 . (25) 

For N independent measurements of the density, we would expect the error 
to decline as iV -1 / 2 . For large enough values of T in dir), where T in d becomes 
independent of r, we can therefore interpret Tind as the time-interval required 
for successive measurements to be independent. 
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The next step of the calculation is to compute the correlation function C(t) 
appropriate for our model. For pure diffusion, we expect: 



C(t) - 1 for t« (26) 

(Ax) 2 \ d/2 , 4 (Ax) 2 



c{t) ~\ m) for t:$> D- (27) 



On time scales t <C (Ax) /D the system remains perfectly correlated as there 
has been insufficient time for particles to hop away to neighboring sites. How- 
ever, for t ^> (Ax) 2 /D, an algebraically decaying correlation function is found, 
characteristic of diffusion. However, we also need to incorporate the effects of 
spontaneous decay that occur independently of the diffusive motion. Adding de- 
cay to the system simply alters the correlation functions by a multiplicative fac- 
tor of exp(-fit). We now substitute this full form into the definition of T% n d ([24)) . 
In the biologically relevant limits where r 3> (Ax) 2 /D and 1/fi (Ax) 2 /D, 
we find, for d = 2 



^( ln ((^) +COnStant ) (28) 



for /ir « 1, and 



r md ~ (in (^a) + constant) (29) 

for /it > 1. In a three-dimensional system we find 

r ind ~^^. (30) 

For the parameter values considered in our simulations we do not observe the 
logarithmic r-dependence in the width predicted by (|28p. In the single gradient 
simulations this is because, at short times r<T x ,we enter the constant w ~ A 
regime. For the parameter values used, the transition from w ~ A at r <C 
t x ss 0.3s to the long time behaviour iJTJl for t » « 4s overwhelms the 
small logarithmic effect. If the production rate J were increased significantly, 
t x oc J -1 would be reduced and the lnr regime would become accessible since 
the r x and l//x timescales would then become better separated. However, even 
in this case, the logarithmic variation in (|28p is intrinsically weak, and will likely 
have a negligible effect in a biological context. 



Simulations 

Stochastic simulations were performed on a two-dimensional square lattice with 
N x = L/Sx sites in the x-direction and N y = Lj_/6x sites in the y-direction, 
where Sx = 0.01/im is the lattice spacing. The detector size Ax was normally 
set equal to 5x except for cases where the detector size was varied, in which 
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case Ax was set to be a multiple of Sx. Zero-flux boundaries were implemented 
at x — and x = L, and a periodic boundary was used to connect y = with 
y = L±. A fixed time step, St = 2.5 x 10 -5 s, was chosen so that for the given 
diffusion constant the total probability of diffusion out of a site in all directions 
approached 1 . However, a timestep 5 times smaller was also tested with no effect 
on any of the results. For each x = site, particles were injected at each time 
step in a Poisson process with mean j = JSxSt. In the two-gradient model, 
particles were also added at x — L in an identical but uncorrelated process. 
Diffusion and decay were also treated as Poisson processes, with hopping and 
decay probabilities of DSt/(5x) 2 and [i8t per particle respectively. Simulations 
were initialised with the mean number of particles in the system, JL±/fi for 
the one-gradient model or twice this value for the two-gradient model, with a 
probability distribution that followed the average density distribution. 

The mean occupancy for each detector site was calculated over the aver- 
aging period, r. For each site this mean occupancy was compared with each 
neighbouring site. If one occupancy was above the threshold and the other be- 
low, this boundary was identified as a threshold crossing position. This process 
was repeated for many averaging periods, ranging from 10 5 repeats for short 
averaging times to 500 repeats for very long averaging times, to generate a dis- 
tribution of crossing positions throughout the system. Threshold crossings in 
both the x- and y— directions were observed. We found that the distributions 
as a function of x— position of these two types of crossing were the same. For 
each row of sites, x = to x — L at fixed y, the mean ("measured threshold") 
and root-mean-squared deviation ( "width" ) of the threshold distribution from 
many averaging periods were calculated independently. In the figures we plot 
the mean of these two quantities across the different y-values within the system, 
with error bars of one standard deviation. 

For the single-source model the standard parameter values used in the sim- 
ulations were as follows: L = 10/im, L± = 6/im, D = l/im 2 s -1 , fi ~ 0.25s -1 , 
J = 4.17/im -1 s -1 , Air = 0.01/im, x^ = 2/im. To generate the data collapse 
in figures [Tp and F, simulations were also performed with: D = 0.5/im 2 s -1 ; 
J = 6.25/im -1 s -1 ; Ax = 0.02/im; fi = Is -1 ; /i = 0.11s -1 ; xt = 1/xm; 
xt = 3/im. For the two-source model, standard parameters were the same as 
above except /i = 0.36s -1 and J = 6/im -1 s -1 . In figure[2p data are also shown 
with: D = 0.5^m 2 s -1 ; ^ = Is -1 ; /i = 0.25s -1 ; J = 9^m -1 s -1 ; Ax = 0.02/^m; 
L = 7.5/im; L = 15/im and Aa: = 0.02/im. 
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Figure 1: Single gradient model in 2c?. A. Variation of the estimated threshold 
position with averaging time, with xt — 2/iin and A = 2/im. B. Variation of the 
width as a function of averaging time. C. Data collapse of the width at large r for 
a range of parameter values. Full line shows the prediction of equation ([7]) with 
k2d = 0.40 and a — 2.5. D. w(t) as a function of decay length, with xt = 2^m. 
Results for three different averaging times are shown: x: t = 10s; o: r = 15s; 
and +: t = 22.5s. The full line shows the prediction from equation ([7]). At large 
A the simulation results deviate from the prediction since the assumption that 
L ^> A is no longer valid. E. Plot of the probability distribution for measuring 
the threshold at position x with an averaging time r = 45s. The full line 
shows a normal distribution. F. Scaling of the cross-over time, r x , according 
to equation (|13j) . In figures A., B. and E. the standard parameter values given 
in the text were used. In figures C. and F., * indicates the standard parameter 
values. For the other data sets one parameter value was changed as follows: o: 
D = O.S/inr^- 1 ; □: J = e^S/xm^s" 1 ; x: Ax = 0.02^m; •: n = Is" 1 ; +: 
H = 0.11s ; o: xt = 1/im; V: xt — 3/im. 
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Figure 2: Two gradient model in 2c?. A. The mean threshold position fluctuates 
about L/2 due to the symmetry of the system. B. Variation of the width w 
as a function of averaging time. C. Data collapse of the width as a function 
of averaging time, at long times, for a range of parameter values. The full line 
shows (|T9|) with k2d = 0.63 and a — 2.5. * indicates the standard parameter 
values. For the other data sets parameter values were changed as follows: o: 
D = Q.b^s- 1 ; □: J = 9jtim _1 a- 1 ; x: Ax = 0.02^m; •: fj, = Is" 1 ; +: 
H = 0.25s- 1 ; o: L = 7.5^m; V: L = 15/im and Ax = 0.02/im. D. Plot of width 
as a function of decay length for averaging times x: r = 30s; o: r = 45s; and 
+ : t = 60s. The full line shows the prediction from equation (|19p. 
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